/*
*  Ambiguity Model   
* -----------------------------------------------------
*
*  Author: Lukas Stoetzer  
*  Date: 26 October 2016
*
*/

data {

  int<lower=1> K; // # items
  int<lower=1> J; // # parties
  int<lower=1> N; // # responses
  int<lower=1> I; //  # experts

  int<lower=0> T; //  # categories

  int<lower=1,upper=K> kk[N]; // # items for n
  int<lower=1,upper=J> jj[N]; // # party for n
  int<lower=1,upper=I> ii[N]; // # expert for n

  int<lower=0> y[N];  // # response for n; y in {1 ... T}

}


parameters {

  ordered[T-1] tau; // # threshold paramters

  vector<lower=0>[K] bk_free; // # free item discrimination
  vector[K-1] ak_free; // # free item  difficulty

  // real<lower=0> sigma_bk; // # variance on bks
  // real<lower=0> sigma_ak; // # variance on aks

  vector[J] x;  // # party platform
  vector<lower=0>[J] nu_free;  // # party ambiguity

}

transformed parameters {

  vector[K] ak;           // constrained item parameters
  vector[K] bk;           // constrained item parameters
  // vector[T-1] tau;           // constrained item parameters
  vector[J] nu;           // constrained item parameters



  // Product Constraints
  bk = bk_free * pow(exp(sum(log(bk_free))), (-inv(K)));
  nu = nu_free * pow(exp(sum(log(nu_free))), (-inv(J)));

  // Sum Constraints
  ak = append_row(ak_free, rep_vector(-1*sum(ak_free), 1));


}

model {
  
    // Auxiallary vectors and reals
    vector[T] prob;  // probabailites
    real eta; // linea r predictor
    real si;   // denominator of linear predicator

    // Prior Specifications
    x ~ normal(0, 1); // Prior on platform with fixed hyper paramters for identification.
    nu_free ~ lognormal(0, 1); // Prior on ambguity with expecattion of 1.

    bk_free ~ lognormal(0, 1 ); // Hierarchical Prior on discrimination; soft identification with fixed mean
    ak_free ~ normal(0, 1); // Hierarchical Prior on difficulty; soft identification with fixed mean

    tau ~ normal(0,1); // prior on threseholds

    // Model; loop over all observations
    for (n in 1:N){

      // Model for latent trait
      eta = bk[kk[n]]* x[jj[n]] + ak[kk[n]];
      si = nu[jj[n]];
      
      // Odered probabilities
      prob[1] = Phi_approx ((tau[1] - eta)/si);
      for (t in 2:(T-1)){
        prob[t] =  Phi_approx((tau[t] - eta)/si) - Phi_approx((tau[t-1] - eta)/si);
      }
      prob[T] =  1- Phi_approx((tau[T-1] - eta)/si);
   
      // Observed outcome
      y[n] ~ categorical(prob);

    }

} 
